In [2]:
import tables
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [3]:
f = tables.open_file("gamma_cone10_sample.h5", "r")

Provenance information stored as attributes

  • Versions of involved code
  • Corsika and simtel main parameters
  • Simtel filenames
In [4]:
f.root._v_attrs
Out[4]:
/._v_attrs (AttributeSet), 20 attributes:
   [CLASS := 'GROUP',
    CORSIKA_ver := None,
    E_max := None,
    E_min := None,
    ImageExtractor_ver := '0.2',
    PYTABLES_FORMAT_VERSION := '2.1',
    TITLE := 'Output File',
    VERSION := '1.0',
    azimuth := 0.0,
    ctapipe_ver := '0.5.2.post312',
    particle_type := 0,
    prod_site_B_field := None,
    prod_site_alt := None,
    prod_site_array := None,
    prod_site_coord := None,
    prod_site_subarray := None,
    runlist := ['/data2/deeplearning/test/gamma_20deg_0deg_run100___cta-prod3-sct_desert-2150m-Paranal-HB9-all_cone10.simtel.gz'],
    simtel_ver := None,
    spectral_index := None,
    zenith := 1.2217304706573486]

Flat file structure: everything else is stored as tables

  • Three tables containing:
    • Array information
    • Telescope type information
    • Event information
  • Plus one table per telescope containing:
    • 2x 1D arrays with calibrated pixels data per image (charge, arrival time)
In [5]:
for node in f:
    print(node)
/ (RootGroup) 'Output File'
/Array_Info (Table(289,)) 'Table of array data'
/Event_Info (Table(100,)) 'Table of Events'
/LST (Table(57,)) 'Table of LST images'
/MSTF (Table(243,)) 'Table of MSTF images'
/MSTN (Table(207,)) 'Table of MSTN images'
/MSTS (Table(188,)) 'Table of MSTS images'
/SST1 (Table(117,)) 'Table of SST1 images'
/SSTA (Table(125,)) 'Table of SSTA images'
/SSTC (Table(123,)) 'Table of SSTC images'
/Telescope_Info (Table(7,)) 'Table of telescope data'

Array information table

In [6]:
f.root.Array_Info
Out[6]:
/Array_Info (Table(289,)) 'Table of array data'
  description := {
  "run_array_direction": Float32Col(shape=(2,), dflt=0.0, pos=0),
  "tel_id": UInt8Col(shape=(), dflt=0, pos=1),
  "tel_type": StringCol(itemsize=8, shape=(), dflt=b'', pos=2),
  "tel_x": Float32Col(shape=(), dflt=0.0, pos=3),
  "tel_y": Float32Col(shape=(), dflt=0.0, pos=4),
  "tel_z": Float32Col(shape=(), dflt=0.0, pos=5)}
  byteorder := 'little'
  chunkshape := (2259,)
In [7]:
tel_types=['LST','MSTS','SSTC']

arr_table=f.root.Array_Info

fig, ax = plt.subplots(figsize=(10,10))

for tel_type in tel_types:
    tel_x = [x['tel_x'] for x in arr_table.iterrows() if x['tel_type'] == tel_type.encode('ascii')]
    tel_y = [x['tel_y'] for x in arr_table.iterrows() if x['tel_type'] == tel_type.encode('ascii')]
    plt.scatter(tel_x, tel_y, label=tel_type)
    
ax.legend()
ax.grid()
plt.show()

Telescope information table

In [8]:
f.root.Telescope_Info
Out[8]:
/Telescope_Info (Table(7,)) 'Table of telescope data'
  description := {
  "num_pixels": UInt32Col(shape=(), dflt=0, pos=0),
  "optical_foclen": Float32Col(shape=(), dflt=0.0, pos=1),
  "tel_type": StringCol(itemsize=8, shape=(), dflt=b'', pos=2),
  "pixel_pos": Float32Col(shape=(2, 11328), dflt=0.0, pos=3)}
  byteorder := 'little'
  chunkshape := (2,)

Event information table

In [9]:
f.root.Event_Info
Out[9]:
/Event_Info (Table(100,)) 'Table of Events'
  description := {
  "alt": Float32Col(shape=(), dflt=0.0, pos=0),
  "az": Float32Col(shape=(), dflt=0.0, pos=1),
  "core_x": Float32Col(shape=(), dflt=0.0, pos=2),
  "core_y": Float32Col(shape=(), dflt=0.0, pos=3),
  "event_number": UInt32Col(shape=(), dflt=0, pos=4),
  "h_first_int": Float32Col(shape=(), dflt=0.0, pos=5),
  "mc_energy": Float32Col(shape=(), dflt=0.0, pos=6),
  "particle_id": UInt8Col(shape=(), dflt=0, pos=7),
  "run_number": UInt32Col(shape=(), dflt=0, pos=8),
  "LST_indices": Int32Col(shape=(4,), dflt=0, pos=9),
  "MSTF_indices": Int32Col(shape=(25,), dflt=0, pos=10),
  "MSTN_indices": Int32Col(shape=(25,), dflt=0, pos=11),
  "MSTS_indices": Int32Col(shape=(25,), dflt=0, pos=12),
  "SST1_indices": Int32Col(shape=(70,), dflt=0, pos=13),
  "SSTA_indices": Int32Col(shape=(70,), dflt=0, pos=14),
  "SSTC_indices": Int32Col(shape=(70,), dflt=0, pos=15)}
  byteorder := 'little'
  chunkshape := (110,)

Example 1: extract some infos from a given event

In [10]:
event_index = 44
my_event = f.root.Event_Info[event_index]
print('Event number: {}'.format(my_event['event_number']))
print('Energy: {} TeV'.format(my_event['mc_energy']))
print('Alt: {} rad'.format(my_event['alt']))
print('Az: {} rad'.format(my_event['az']))
tel_types = ['LST','MSTS','SSTC']
for tel_type in tel_types:
    tel_indices = '{}_indices'.format(tel_type)
    print('{} = {}'.format(tel_indices,my_event[tel_indices]))
Event number: 296407
Energy: 0.6381028294563293 TeV
Alt: 1.2547309398651123 rad
Az: 0.0808204710483551 rad
LST_indices = [ 0 22  0  0]
MSTS_indices = [76  0 77 78 79  0 80  0 81  0  0  0  0  0 82  0 83  0  0  0  0  0  0  0  0]
SSTC_indices = [ 0 53  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]

Example 2: plot all charge and arrival time images for a given event

In [11]:
tel_types=['LST','MSTS','MSTN','MSTF','SSTC','SST1','SSTA']

tel_type_size={'LST':15,'MSTS':2,'MSTN':15,'MSTF':15,'SSTC':15,'SSTA':15,'SST1':15} #Size of the scatter plot point, for visualization purposes only

tel_table=f.root.Telescope_Info

for tel_type in tel_types:
    pos_x, pos_y = [x['pixel_pos'] for x in tel_table.iterrows() if x['tel_type'] == tel_type.encode('ascii')][0]
    tel_ids = [x['tel_id'] for x in arr_table.iterrows() if x['tel_type'] == tel_type.encode('ascii')]
    exec('my_images = f.root.{}'.format(tel_type))
    my_indices = my_event['{}_indices'.format(tel_type)]
    for img_index in my_indices:
        if img_index > 0:
            img_charge = my_images[img_index]['image_charge']
            img_time = my_images[img_index]['image_peak_times']
            print('Telescope type: {}, ID: {}'.format(tel_type,tel_ids[list(my_indices).index(img_index)]))
            plt.figure(figsize=(15,5))
            plt.subplot(1, 2, 1)
            plt.scatter(pos_x[:len(img_charge)], pos_y[:len(img_charge)], c=img_charge, s=tel_type_size[tel_type])
            plt.axis('square')
            plt.colorbar()
            plt.subplot(1, 2, 2)
            plt.scatter(pos_x[:len(img_time)], pos_y[:len(img_time)], c=img_time, s=tel_type_size[tel_type])
            plt.axis('square')
            plt.colorbar()
            plt.show()
Telescope type: LST, ID: 2
Telescope type: MSTS, ID: 55
Telescope type: MSTS, ID: 57
Telescope type: MSTS, ID: 58
Telescope type: MSTS, ID: 59
Telescope type: MSTS, ID: 61
Telescope type: MSTS, ID: 63
Telescope type: MSTS, ID: 69
Telescope type: MSTS, ID: 71
Telescope type: MSTN, ID: 5
Telescope type: MSTN, ID: 7
Telescope type: MSTN, ID: 8
Telescope type: MSTN, ID: 9
Telescope type: MSTN, ID: 11
Telescope type: MSTN, ID: 13
Telescope type: MSTN, ID: 16
Telescope type: MSTN, ID: 19
Telescope type: MSTN, ID: 21
Telescope type: MSTF, ID: 30
Telescope type: MSTF, ID: 32
Telescope type: MSTF, ID: 34
Telescope type: MSTF, ID: 36
Telescope type: MSTF, ID: 38
Telescope type: MSTF, ID: 41
Telescope type: MSTF, ID: 44
Telescope type: MSTF, ID: 46
Telescope type: MSTF, ID: 50
Telescope type: SSTC, ID: 151
Telescope type: SST1, ID: 221
Telescope type: SSTA, ID: 81